Diagnostics and Posteriors

Elizabeth King
Kevin Middleton

Learning path

rethinking package (McElreath 2020)

  • Between pure stan code and convenience tools
  • Teaches careful practices, diagnosis, summaries
  • Lacks integration with other tools

“Automated” packages: rstanarm and brms

  • Work well with other tools (posterior, tidybayes, bayesplot)
  • Graduate to these packages once you are comfortable

Diagnostics

  • Bayesian inference is active
    • More than just lm() and summary()
  • Define a model
  • Prior prediction
  • Sampling
  • Check that sampling was “good”
  • Posterior prediction

Use our simulated data

set.seed(467456)
(y <- rnorm(n = 8, mean = 50, sd = 2))
[1] 44.88692 48.25894 48.44770 50.16844 49.55441 49.64544 48.49286 48.91631
mean(y)
[1] 48.54638
sd(y)
[1] 1.624759

Write the model statements

Likelihood:

\[y \sim Normal(\mu, \sigma)\]

Priors:

\[\mu \sim Normal(45, 10)\]

\[\sigma \sim HalfNormal(0, 3)\]

Look at the priors

Skip prior predictive simulation (for now)

Create ulam() model

Note the correspondence to the model statements.

fm <- ulam(
  alist(
    y ~ dnorm(mu, sigma),
    mu ~ dnorm(45, 10),
    sigma ~ dhalfnorm(0, 3)
  ),
  data = list(y = y),
  chains = 4,
  iter = 1e4
)

50% of iterations are used for warm-up (default)

Create ulam() model

Running MCMC with 4 sequential chains, with 1 thread(s) per chain...

Chain 1 Iteration:    1 / 10000 [  0%]  (Warmup) 
Chain 1 Iteration:  100 / 10000 [  1%]  (Warmup) 
Chain 1 Iteration:  200 / 10000 [  2%]  (Warmup) 
Chain 1 Iteration:  300 / 10000 [  3%]  (Warmup) 
Chain 1 Iteration:  400 / 10000 [  4%]  (Warmup) 
Chain 1 Iteration:  500 / 10000 [  5%]  (Warmup) 
Chain 1 Iteration:  600 / 10000 [  6%]  (Warmup) 
Chain 1 Iteration:  700 / 10000 [  7%]  (Warmup) 
Chain 1 Iteration:  800 / 10000 [  8%]  (Warmup) 
Chain 1 Iteration:  900 / 10000 [  9%]  (Warmup) 
Chain 1 Iteration: 1000 / 10000 [ 10%]  (Warmup) 
Chain 1 Iteration: 1100 / 10000 [ 11%]  (Warmup) 
Chain 1 Iteration: 1200 / 10000 [ 12%]  (Warmup) 
Chain 1 Iteration: 1300 / 10000 [ 13%]  (Warmup) 
Chain 1 Iteration: 1400 / 10000 [ 14%]  (Warmup) 
Chain 1 Iteration: 1500 / 10000 [ 15%]  (Warmup) 
Chain 1 Iteration: 1600 / 10000 [ 16%]  (Warmup) 
Chain 1 Iteration: 1700 / 10000 [ 17%]  (Warmup) 
Chain 1 Iteration: 1800 / 10000 [ 18%]  (Warmup) 
Chain 1 Iteration: 1900 / 10000 [ 19%]  (Warmup) 
Chain 1 Iteration: 2000 / 10000 [ 20%]  (Warmup) 
Chain 1 Iteration: 2100 / 10000 [ 21%]  (Warmup) 
Chain 1 Iteration: 2200 / 10000 [ 22%]  (Warmup) 
Chain 1 Iteration: 2300 / 10000 [ 23%]  (Warmup) 
Chain 1 Iteration: 2400 / 10000 [ 24%]  (Warmup) 
Chain 1 Iteration: 2500 / 10000 [ 25%]  (Warmup) 
Chain 1 Iteration: 2600 / 10000 [ 26%]  (Warmup) 
Chain 1 Iteration: 2700 / 10000 [ 27%]  (Warmup) 
Chain 1 Iteration: 2800 / 10000 [ 28%]  (Warmup) 
Chain 1 Iteration: 2900 / 10000 [ 29%]  (Warmup) 
Chain 1 Iteration: 3000 / 10000 [ 30%]  (Warmup) 
Chain 1 Iteration: 3100 / 10000 [ 31%]  (Warmup) 
Chain 1 Iteration: 3200 / 10000 [ 32%]  (Warmup) 
Chain 1 Iteration: 3300 / 10000 [ 33%]  (Warmup) 
Chain 1 Iteration: 3400 / 10000 [ 34%]  (Warmup) 
Chain 1 Iteration: 3500 / 10000 [ 35%]  (Warmup) 
Chain 1 Iteration: 3600 / 10000 [ 36%]  (Warmup) 
Chain 1 Iteration: 3700 / 10000 [ 37%]  (Warmup) 
Chain 1 Iteration: 3800 / 10000 [ 38%]  (Warmup) 
Chain 1 Iteration: 3900 / 10000 [ 39%]  (Warmup) 
Chain 1 Iteration: 4000 / 10000 [ 40%]  (Warmup) 
Chain 1 Iteration: 4100 / 10000 [ 41%]  (Warmup) 
Chain 1 Iteration: 4200 / 10000 [ 42%]  (Warmup) 
Chain 1 Iteration: 4300 / 10000 [ 43%]  (Warmup) 
Chain 1 Iteration: 4400 / 10000 [ 44%]  (Warmup) 
Chain 1 Iteration: 4500 / 10000 [ 45%]  (Warmup) 
Chain 1 Iteration: 4600 / 10000 [ 46%]  (Warmup) 
Chain 1 Iteration: 4700 / 10000 [ 47%]  (Warmup) 
Chain 1 Iteration: 4800 / 10000 [ 48%]  (Warmup) 
Chain 1 Iteration: 4900 / 10000 [ 49%]  (Warmup) 
Chain 1 Iteration: 5000 / 10000 [ 50%]  (Warmup) 
Chain 1 Iteration: 5001 / 10000 [ 50%]  (Sampling) 
Chain 1 Iteration: 5100 / 10000 [ 51%]  (Sampling) 
Chain 1 Iteration: 5200 / 10000 [ 52%]  (Sampling) 
Chain 1 Iteration: 5300 / 10000 [ 53%]  (Sampling) 
Chain 1 Iteration: 5400 / 10000 [ 54%]  (Sampling) 
Chain 1 Iteration: 5500 / 10000 [ 55%]  (Sampling) 
Chain 1 Iteration: 5600 / 10000 [ 56%]  (Sampling) 
Chain 1 Iteration: 5700 / 10000 [ 57%]  (Sampling) 
Chain 1 Iteration: 5800 / 10000 [ 58%]  (Sampling) 
Chain 1 Iteration: 5900 / 10000 [ 59%]  (Sampling) 
Chain 1 Iteration: 6000 / 10000 [ 60%]  (Sampling) 
Chain 1 Iteration: 6100 / 10000 [ 61%]  (Sampling) 
Chain 1 Iteration: 6200 / 10000 [ 62%]  (Sampling) 
Chain 1 Iteration: 6300 / 10000 [ 63%]  (Sampling) 
Chain 1 Iteration: 6400 / 10000 [ 64%]  (Sampling) 
Chain 1 Iteration: 6500 / 10000 [ 65%]  (Sampling) 
Chain 1 Iteration: 6600 / 10000 [ 66%]  (Sampling) 
Chain 1 Iteration: 6700 / 10000 [ 67%]  (Sampling) 
Chain 1 Iteration: 6800 / 10000 [ 68%]  (Sampling) 
Chain 1 Iteration: 6900 / 10000 [ 69%]  (Sampling) 
Chain 1 Iteration: 7000 / 10000 [ 70%]  (Sampling) 
Chain 1 Iteration: 7100 / 10000 [ 71%]  (Sampling) 
Chain 1 Iteration: 7200 / 10000 [ 72%]  (Sampling) 
Chain 1 Iteration: 7300 / 10000 [ 73%]  (Sampling) 
Chain 1 Iteration: 7400 / 10000 [ 74%]  (Sampling) 
Chain 1 Iteration: 7500 / 10000 [ 75%]  (Sampling) 
Chain 1 Iteration: 7600 / 10000 [ 76%]  (Sampling) 
Chain 1 Iteration: 7700 / 10000 [ 77%]  (Sampling) 
Chain 1 Iteration: 7800 / 10000 [ 78%]  (Sampling) 
Chain 1 Iteration: 7900 / 10000 [ 79%]  (Sampling) 
Chain 1 Iteration: 8000 / 10000 [ 80%]  (Sampling) 
Chain 1 Iteration: 8100 / 10000 [ 81%]  (Sampling) 
Chain 1 Iteration: 8200 / 10000 [ 82%]  (Sampling) 
Chain 1 Iteration: 8300 / 10000 [ 83%]  (Sampling) 
Chain 1 Iteration: 8400 / 10000 [ 84%]  (Sampling) 
Chain 1 Iteration: 8500 / 10000 [ 85%]  (Sampling) 
Chain 1 Iteration: 8600 / 10000 [ 86%]  (Sampling) 
Chain 1 Iteration: 8700 / 10000 [ 87%]  (Sampling) 
Chain 1 Iteration: 8800 / 10000 [ 88%]  (Sampling) 
Chain 1 Iteration: 8900 / 10000 [ 89%]  (Sampling) 
Chain 1 Iteration: 9000 / 10000 [ 90%]  (Sampling) 
Chain 1 Iteration: 9100 / 10000 [ 91%]  (Sampling) 
Chain 1 Iteration: 9200 / 10000 [ 92%]  (Sampling) 
Chain 1 Iteration: 9300 / 10000 [ 93%]  (Sampling) 
Chain 1 Iteration: 9400 / 10000 [ 94%]  (Sampling) 
Chain 1 Iteration: 9500 / 10000 [ 95%]  (Sampling) 
Chain 1 Iteration: 9600 / 10000 [ 96%]  (Sampling) 
Chain 1 Iteration: 9700 / 10000 [ 97%]  (Sampling) 
Chain 1 Iteration: 9800 / 10000 [ 98%]  (Sampling) 
Chain 1 Iteration: 9900 / 10000 [ 99%]  (Sampling) 
Chain 1 Iteration: 10000 / 10000 [100%]  (Sampling) 
Chain 1 finished in 0.1 seconds.
Chain 2 Iteration:    1 / 10000 [  0%]  (Warmup) 
Chain 2 Iteration:  100 / 10000 [  1%]  (Warmup) 
Chain 2 Iteration:  200 / 10000 [  2%]  (Warmup) 
Chain 2 Iteration:  300 / 10000 [  3%]  (Warmup) 
Chain 2 Iteration:  400 / 10000 [  4%]  (Warmup) 
Chain 2 Iteration:  500 / 10000 [  5%]  (Warmup) 
Chain 2 Iteration:  600 / 10000 [  6%]  (Warmup) 
Chain 2 Iteration:  700 / 10000 [  7%]  (Warmup) 
Chain 2 Iteration:  800 / 10000 [  8%]  (Warmup) 
Chain 2 Iteration:  900 / 10000 [  9%]  (Warmup) 
Chain 2 Iteration: 1000 / 10000 [ 10%]  (Warmup) 
Chain 2 Iteration: 1100 / 10000 [ 11%]  (Warmup) 
Chain 2 Iteration: 1200 / 10000 [ 12%]  (Warmup) 
Chain 2 Iteration: 1300 / 10000 [ 13%]  (Warmup) 
Chain 2 Iteration: 1400 / 10000 [ 14%]  (Warmup) 
Chain 2 Iteration: 1500 / 10000 [ 15%]  (Warmup) 
Chain 2 Iteration: 1600 / 10000 [ 16%]  (Warmup) 
Chain 2 Iteration: 1700 / 10000 [ 17%]  (Warmup) 
Chain 2 Iteration: 1800 / 10000 [ 18%]  (Warmup) 
Chain 2 Iteration: 1900 / 10000 [ 19%]  (Warmup) 
Chain 2 Iteration: 2000 / 10000 [ 20%]  (Warmup) 
Chain 2 Iteration: 2100 / 10000 [ 21%]  (Warmup) 
Chain 2 Iteration: 2200 / 10000 [ 22%]  (Warmup) 
Chain 2 Iteration: 2300 / 10000 [ 23%]  (Warmup) 
Chain 2 Iteration: 2400 / 10000 [ 24%]  (Warmup) 
Chain 2 Iteration: 2500 / 10000 [ 25%]  (Warmup) 
Chain 2 Iteration: 2600 / 10000 [ 26%]  (Warmup) 
Chain 2 Iteration: 2700 / 10000 [ 27%]  (Warmup) 
Chain 2 Iteration: 2800 / 10000 [ 28%]  (Warmup) 
Chain 2 Iteration: 2900 / 10000 [ 29%]  (Warmup) 
Chain 2 Iteration: 3000 / 10000 [ 30%]  (Warmup) 
Chain 2 Iteration: 3100 / 10000 [ 31%]  (Warmup) 
Chain 2 Iteration: 3200 / 10000 [ 32%]  (Warmup) 
Chain 2 Iteration: 3300 / 10000 [ 33%]  (Warmup) 
Chain 2 Iteration: 3400 / 10000 [ 34%]  (Warmup) 
Chain 2 Iteration: 3500 / 10000 [ 35%]  (Warmup) 
Chain 2 Iteration: 3600 / 10000 [ 36%]  (Warmup) 
Chain 2 Iteration: 3700 / 10000 [ 37%]  (Warmup) 
Chain 2 Iteration: 3800 / 10000 [ 38%]  (Warmup) 
Chain 2 Iteration: 3900 / 10000 [ 39%]  (Warmup) 
Chain 2 Iteration: 4000 / 10000 [ 40%]  (Warmup) 
Chain 2 Iteration: 4100 / 10000 [ 41%]  (Warmup) 
Chain 2 Iteration: 4200 / 10000 [ 42%]  (Warmup) 
Chain 2 Iteration: 4300 / 10000 [ 43%]  (Warmup) 
Chain 2 Iteration: 4400 / 10000 [ 44%]  (Warmup) 
Chain 2 Iteration: 4500 / 10000 [ 45%]  (Warmup) 
Chain 2 Iteration: 4600 / 10000 [ 46%]  (Warmup) 
Chain 2 Iteration: 4700 / 10000 [ 47%]  (Warmup) 
Chain 2 Iteration: 4800 / 10000 [ 48%]  (Warmup) 
Chain 2 Iteration: 4900 / 10000 [ 49%]  (Warmup) 
Chain 2 Iteration: 5000 / 10000 [ 50%]  (Warmup) 
Chain 2 Iteration: 5001 / 10000 [ 50%]  (Sampling) 
Chain 2 Iteration: 5100 / 10000 [ 51%]  (Sampling) 
Chain 2 Iteration: 5200 / 10000 [ 52%]  (Sampling) 
Chain 2 Iteration: 5300 / 10000 [ 53%]  (Sampling) 
Chain 2 Iteration: 5400 / 10000 [ 54%]  (Sampling) 
Chain 2 Iteration: 5500 / 10000 [ 55%]  (Sampling) 
Chain 2 Iteration: 5600 / 10000 [ 56%]  (Sampling) 
Chain 2 Iteration: 5700 / 10000 [ 57%]  (Sampling) 
Chain 2 Iteration: 5800 / 10000 [ 58%]  (Sampling) 
Chain 2 Iteration: 5900 / 10000 [ 59%]  (Sampling) 
Chain 2 Iteration: 6000 / 10000 [ 60%]  (Sampling) 
Chain 2 Iteration: 6100 / 10000 [ 61%]  (Sampling) 
Chain 2 Iteration: 6200 / 10000 [ 62%]  (Sampling) 
Chain 2 Iteration: 6300 / 10000 [ 63%]  (Sampling) 
Chain 2 Iteration: 6400 / 10000 [ 64%]  (Sampling) 
Chain 2 Iteration: 6500 / 10000 [ 65%]  (Sampling) 
Chain 2 Iteration: 6600 / 10000 [ 66%]  (Sampling) 
Chain 2 Iteration: 6700 / 10000 [ 67%]  (Sampling) 
Chain 2 Iteration: 6800 / 10000 [ 68%]  (Sampling) 
Chain 2 Iteration: 6900 / 10000 [ 69%]  (Sampling) 
Chain 2 Iteration: 7000 / 10000 [ 70%]  (Sampling) 
Chain 2 Iteration: 7100 / 10000 [ 71%]  (Sampling) 
Chain 2 Iteration: 7200 / 10000 [ 72%]  (Sampling) 
Chain 2 Iteration: 7300 / 10000 [ 73%]  (Sampling) 
Chain 2 Iteration: 7400 / 10000 [ 74%]  (Sampling) 
Chain 2 Iteration: 7500 / 10000 [ 75%]  (Sampling) 
Chain 2 Iteration: 7600 / 10000 [ 76%]  (Sampling) 
Chain 2 Iteration: 7700 / 10000 [ 77%]  (Sampling) 
Chain 2 Iteration: 7800 / 10000 [ 78%]  (Sampling) 
Chain 2 Iteration: 7900 / 10000 [ 79%]  (Sampling) 
Chain 2 Iteration: 8000 / 10000 [ 80%]  (Sampling) 
Chain 2 Iteration: 8100 / 10000 [ 81%]  (Sampling) 
Chain 2 Iteration: 8200 / 10000 [ 82%]  (Sampling) 
Chain 2 Iteration: 8300 / 10000 [ 83%]  (Sampling) 
Chain 2 Iteration: 8400 / 10000 [ 84%]  (Sampling) 
Chain 2 Iteration: 8500 / 10000 [ 85%]  (Sampling) 
Chain 2 Iteration: 8600 / 10000 [ 86%]  (Sampling) 
Chain 2 Iteration: 8700 / 10000 [ 87%]  (Sampling) 
Chain 2 Iteration: 8800 / 10000 [ 88%]  (Sampling) 
Chain 2 Iteration: 8900 / 10000 [ 89%]  (Sampling) 
Chain 2 Iteration: 9000 / 10000 [ 90%]  (Sampling) 
Chain 2 Iteration: 9100 / 10000 [ 91%]  (Sampling) 
Chain 2 Iteration: 9200 / 10000 [ 92%]  (Sampling) 
Chain 2 Iteration: 9300 / 10000 [ 93%]  (Sampling) 
Chain 2 Iteration: 9400 / 10000 [ 94%]  (Sampling) 
Chain 2 Iteration: 9500 / 10000 [ 95%]  (Sampling) 
Chain 2 Iteration: 9600 / 10000 [ 96%]  (Sampling) 
Chain 2 Iteration: 9700 / 10000 [ 97%]  (Sampling) 
Chain 2 Iteration: 9800 / 10000 [ 98%]  (Sampling) 
Chain 2 Iteration: 9900 / 10000 [ 99%]  (Sampling) 
Chain 2 Iteration: 10000 / 10000 [100%]  (Sampling) 
Chain 2 finished in 0.1 seconds.
Chain 3 Iteration:    1 / 10000 [  0%]  (Warmup) 
Chain 3 Iteration:  100 / 10000 [  1%]  (Warmup) 
Chain 3 Iteration:  200 / 10000 [  2%]  (Warmup) 
Chain 3 Iteration:  300 / 10000 [  3%]  (Warmup) 
Chain 3 Iteration:  400 / 10000 [  4%]  (Warmup) 
Chain 3 Iteration:  500 / 10000 [  5%]  (Warmup) 
Chain 3 Iteration:  600 / 10000 [  6%]  (Warmup) 
Chain 3 Iteration:  700 / 10000 [  7%]  (Warmup) 
Chain 3 Iteration:  800 / 10000 [  8%]  (Warmup) 
Chain 3 Iteration:  900 / 10000 [  9%]  (Warmup) 
Chain 3 Iteration: 1000 / 10000 [ 10%]  (Warmup) 
Chain 3 Iteration: 1100 / 10000 [ 11%]  (Warmup) 
Chain 3 Iteration: 1200 / 10000 [ 12%]  (Warmup) 
Chain 3 Iteration: 1300 / 10000 [ 13%]  (Warmup) 
Chain 3 Iteration: 1400 / 10000 [ 14%]  (Warmup) 
Chain 3 Iteration: 1500 / 10000 [ 15%]  (Warmup) 
Chain 3 Iteration: 1600 / 10000 [ 16%]  (Warmup) 
Chain 3 Iteration: 1700 / 10000 [ 17%]  (Warmup) 
Chain 3 Iteration: 1800 / 10000 [ 18%]  (Warmup) 
Chain 3 Iteration: 1900 / 10000 [ 19%]  (Warmup) 
Chain 3 Iteration: 2000 / 10000 [ 20%]  (Warmup) 
Chain 3 Iteration: 2100 / 10000 [ 21%]  (Warmup) 
Chain 3 Iteration: 2200 / 10000 [ 22%]  (Warmup) 
Chain 3 Iteration: 2300 / 10000 [ 23%]  (Warmup) 
Chain 3 Iteration: 2400 / 10000 [ 24%]  (Warmup) 
Chain 3 Iteration: 2500 / 10000 [ 25%]  (Warmup) 
Chain 3 Iteration: 2600 / 10000 [ 26%]  (Warmup) 
Chain 3 Iteration: 2700 / 10000 [ 27%]  (Warmup) 
Chain 3 Iteration: 2800 / 10000 [ 28%]  (Warmup) 
Chain 3 Iteration: 2900 / 10000 [ 29%]  (Warmup) 
Chain 3 Iteration: 3000 / 10000 [ 30%]  (Warmup) 
Chain 3 Iteration: 3100 / 10000 [ 31%]  (Warmup) 
Chain 3 Iteration: 3200 / 10000 [ 32%]  (Warmup) 
Chain 3 Iteration: 3300 / 10000 [ 33%]  (Warmup) 
Chain 3 Iteration: 3400 / 10000 [ 34%]  (Warmup) 
Chain 3 Iteration: 3500 / 10000 [ 35%]  (Warmup) 
Chain 3 Iteration: 3600 / 10000 [ 36%]  (Warmup) 
Chain 3 Iteration: 3700 / 10000 [ 37%]  (Warmup) 
Chain 3 Iteration: 3800 / 10000 [ 38%]  (Warmup) 
Chain 3 Iteration: 3900 / 10000 [ 39%]  (Warmup) 
Chain 3 Iteration: 4000 / 10000 [ 40%]  (Warmup) 
Chain 3 Iteration: 4100 / 10000 [ 41%]  (Warmup) 
Chain 3 Iteration: 4200 / 10000 [ 42%]  (Warmup) 
Chain 3 Iteration: 4300 / 10000 [ 43%]  (Warmup) 
Chain 3 Iteration: 4400 / 10000 [ 44%]  (Warmup) 
Chain 3 Iteration: 4500 / 10000 [ 45%]  (Warmup) 
Chain 3 Iteration: 4600 / 10000 [ 46%]  (Warmup) 
Chain 3 Iteration: 4700 / 10000 [ 47%]  (Warmup) 
Chain 3 Iteration: 4800 / 10000 [ 48%]  (Warmup) 
Chain 3 Iteration: 4900 / 10000 [ 49%]  (Warmup) 
Chain 3 Iteration: 5000 / 10000 [ 50%]  (Warmup) 
Chain 3 Iteration: 5001 / 10000 [ 50%]  (Sampling) 
Chain 3 Iteration: 5100 / 10000 [ 51%]  (Sampling) 
Chain 3 Iteration: 5200 / 10000 [ 52%]  (Sampling) 
Chain 3 Iteration: 5300 / 10000 [ 53%]  (Sampling) 
Chain 3 Iteration: 5400 / 10000 [ 54%]  (Sampling) 
Chain 3 Iteration: 5500 / 10000 [ 55%]  (Sampling) 
Chain 3 Iteration: 5600 / 10000 [ 56%]  (Sampling) 
Chain 3 Iteration: 5700 / 10000 [ 57%]  (Sampling) 
Chain 3 Iteration: 5800 / 10000 [ 58%]  (Sampling) 
Chain 3 Iteration: 5900 / 10000 [ 59%]  (Sampling) 
Chain 3 Iteration: 6000 / 10000 [ 60%]  (Sampling) 
Chain 3 Iteration: 6100 / 10000 [ 61%]  (Sampling) 
Chain 3 Iteration: 6200 / 10000 [ 62%]  (Sampling) 
Chain 3 Iteration: 6300 / 10000 [ 63%]  (Sampling) 
Chain 3 Iteration: 6400 / 10000 [ 64%]  (Sampling) 
Chain 3 Iteration: 6500 / 10000 [ 65%]  (Sampling) 
Chain 3 Iteration: 6600 / 10000 [ 66%]  (Sampling) 
Chain 3 Iteration: 6700 / 10000 [ 67%]  (Sampling) 
Chain 3 Iteration: 6800 / 10000 [ 68%]  (Sampling) 
Chain 3 Iteration: 6900 / 10000 [ 69%]  (Sampling) 
Chain 3 Iteration: 7000 / 10000 [ 70%]  (Sampling) 
Chain 3 Iteration: 7100 / 10000 [ 71%]  (Sampling) 
Chain 3 Iteration: 7200 / 10000 [ 72%]  (Sampling) 
Chain 3 Iteration: 7300 / 10000 [ 73%]  (Sampling) 
Chain 3 Iteration: 7400 / 10000 [ 74%]  (Sampling) 
Chain 3 Iteration: 7500 / 10000 [ 75%]  (Sampling) 
Chain 3 Iteration: 7600 / 10000 [ 76%]  (Sampling) 
Chain 3 Iteration: 7700 / 10000 [ 77%]  (Sampling) 
Chain 3 Iteration: 7800 / 10000 [ 78%]  (Sampling) 
Chain 3 Iteration: 7900 / 10000 [ 79%]  (Sampling) 
Chain 3 Iteration: 8000 / 10000 [ 80%]  (Sampling) 
Chain 3 Iteration: 8100 / 10000 [ 81%]  (Sampling) 
Chain 3 Iteration: 8200 / 10000 [ 82%]  (Sampling) 
Chain 3 Iteration: 8300 / 10000 [ 83%]  (Sampling) 
Chain 3 Iteration: 8400 / 10000 [ 84%]  (Sampling) 
Chain 3 Iteration: 8500 / 10000 [ 85%]  (Sampling) 
Chain 3 Iteration: 8600 / 10000 [ 86%]  (Sampling) 
Chain 3 Iteration: 8700 / 10000 [ 87%]  (Sampling) 
Chain 3 Iteration: 8800 / 10000 [ 88%]  (Sampling) 
Chain 3 Iteration: 8900 / 10000 [ 89%]  (Sampling) 
Chain 3 Iteration: 9000 / 10000 [ 90%]  (Sampling) 
Chain 3 Iteration: 9100 / 10000 [ 91%]  (Sampling) 
Chain 3 Iteration: 9200 / 10000 [ 92%]  (Sampling) 
Chain 3 Iteration: 9300 / 10000 [ 93%]  (Sampling) 
Chain 3 Iteration: 9400 / 10000 [ 94%]  (Sampling) 
Chain 3 Iteration: 9500 / 10000 [ 95%]  (Sampling) 
Chain 3 Iteration: 9600 / 10000 [ 96%]  (Sampling) 
Chain 3 Iteration: 9700 / 10000 [ 97%]  (Sampling) 
Chain 3 Iteration: 9800 / 10000 [ 98%]  (Sampling) 
Chain 3 Iteration: 9900 / 10000 [ 99%]  (Sampling) 
Chain 3 Iteration: 10000 / 10000 [100%]  (Sampling) 
Chain 3 finished in 0.1 seconds.
Chain 4 Iteration:    1 / 10000 [  0%]  (Warmup) 
Chain 4 Iteration:  100 / 10000 [  1%]  (Warmup) 
Chain 4 Iteration:  200 / 10000 [  2%]  (Warmup) 
Chain 4 Iteration:  300 / 10000 [  3%]  (Warmup) 
Chain 4 Iteration:  400 / 10000 [  4%]  (Warmup) 
Chain 4 Iteration:  500 / 10000 [  5%]  (Warmup) 
Chain 4 Iteration:  600 / 10000 [  6%]  (Warmup) 
Chain 4 Iteration:  700 / 10000 [  7%]  (Warmup) 
Chain 4 Iteration:  800 / 10000 [  8%]  (Warmup) 
Chain 4 Iteration:  900 / 10000 [  9%]  (Warmup) 
Chain 4 Iteration: 1000 / 10000 [ 10%]  (Warmup) 
Chain 4 Iteration: 1100 / 10000 [ 11%]  (Warmup) 
Chain 4 Iteration: 1200 / 10000 [ 12%]  (Warmup) 
Chain 4 Iteration: 1300 / 10000 [ 13%]  (Warmup) 
Chain 4 Iteration: 1400 / 10000 [ 14%]  (Warmup) 
Chain 4 Iteration: 1500 / 10000 [ 15%]  (Warmup) 
Chain 4 Iteration: 1600 / 10000 [ 16%]  (Warmup) 
Chain 4 Iteration: 1700 / 10000 [ 17%]  (Warmup) 
Chain 4 Iteration: 1800 / 10000 [ 18%]  (Warmup) 
Chain 4 Iteration: 1900 / 10000 [ 19%]  (Warmup) 
Chain 4 Iteration: 2000 / 10000 [ 20%]  (Warmup) 
Chain 4 Iteration: 2100 / 10000 [ 21%]  (Warmup) 
Chain 4 Iteration: 2200 / 10000 [ 22%]  (Warmup) 
Chain 4 Iteration: 2300 / 10000 [ 23%]  (Warmup) 
Chain 4 Iteration: 2400 / 10000 [ 24%]  (Warmup) 
Chain 4 Iteration: 2500 / 10000 [ 25%]  (Warmup) 
Chain 4 Iteration: 2600 / 10000 [ 26%]  (Warmup) 
Chain 4 Iteration: 2700 / 10000 [ 27%]  (Warmup) 
Chain 4 Iteration: 2800 / 10000 [ 28%]  (Warmup) 
Chain 4 Iteration: 2900 / 10000 [ 29%]  (Warmup) 
Chain 4 Iteration: 3000 / 10000 [ 30%]  (Warmup) 
Chain 4 Iteration: 3100 / 10000 [ 31%]  (Warmup) 
Chain 4 Iteration: 3200 / 10000 [ 32%]  (Warmup) 
Chain 4 Iteration: 3300 / 10000 [ 33%]  (Warmup) 
Chain 4 Iteration: 3400 / 10000 [ 34%]  (Warmup) 
Chain 4 Iteration: 3500 / 10000 [ 35%]  (Warmup) 
Chain 4 Iteration: 3600 / 10000 [ 36%]  (Warmup) 
Chain 4 Iteration: 3700 / 10000 [ 37%]  (Warmup) 
Chain 4 Iteration: 3800 / 10000 [ 38%]  (Warmup) 
Chain 4 Iteration: 3900 / 10000 [ 39%]  (Warmup) 
Chain 4 Iteration: 4000 / 10000 [ 40%]  (Warmup) 
Chain 4 Iteration: 4100 / 10000 [ 41%]  (Warmup) 
Chain 4 Iteration: 4200 / 10000 [ 42%]  (Warmup) 
Chain 4 Iteration: 4300 / 10000 [ 43%]  (Warmup) 
Chain 4 Iteration: 4400 / 10000 [ 44%]  (Warmup) 
Chain 4 Iteration: 4500 / 10000 [ 45%]  (Warmup) 
Chain 4 Iteration: 4600 / 10000 [ 46%]  (Warmup) 
Chain 4 Iteration: 4700 / 10000 [ 47%]  (Warmup) 
Chain 4 Iteration: 4800 / 10000 [ 48%]  (Warmup) 
Chain 4 Iteration: 4900 / 10000 [ 49%]  (Warmup) 
Chain 4 Iteration: 5000 / 10000 [ 50%]  (Warmup) 
Chain 4 Iteration: 5001 / 10000 [ 50%]  (Sampling) 
Chain 4 Iteration: 5100 / 10000 [ 51%]  (Sampling) 
Chain 4 Iteration: 5200 / 10000 [ 52%]  (Sampling) 
Chain 4 Iteration: 5300 / 10000 [ 53%]  (Sampling) 
Chain 4 Iteration: 5400 / 10000 [ 54%]  (Sampling) 
Chain 4 Iteration: 5500 / 10000 [ 55%]  (Sampling) 
Chain 4 Iteration: 5600 / 10000 [ 56%]  (Sampling) 
Chain 4 Iteration: 5700 / 10000 [ 57%]  (Sampling) 
Chain 4 Iteration: 5800 / 10000 [ 58%]  (Sampling) 
Chain 4 Iteration: 5900 / 10000 [ 59%]  (Sampling) 
Chain 4 Iteration: 6000 / 10000 [ 60%]  (Sampling) 
Chain 4 Iteration: 6100 / 10000 [ 61%]  (Sampling) 
Chain 4 Iteration: 6200 / 10000 [ 62%]  (Sampling) 
Chain 4 Iteration: 6300 / 10000 [ 63%]  (Sampling) 
Chain 4 Iteration: 6400 / 10000 [ 64%]  (Sampling) 
Chain 4 Iteration: 6500 / 10000 [ 65%]  (Sampling) 
Chain 4 Iteration: 6600 / 10000 [ 66%]  (Sampling) 
Chain 4 Iteration: 6700 / 10000 [ 67%]  (Sampling) 
Chain 4 Iteration: 6800 / 10000 [ 68%]  (Sampling) 
Chain 4 Iteration: 6900 / 10000 [ 69%]  (Sampling) 
Chain 4 Iteration: 7000 / 10000 [ 70%]  (Sampling) 
Chain 4 Iteration: 7100 / 10000 [ 71%]  (Sampling) 
Chain 4 Iteration: 7200 / 10000 [ 72%]  (Sampling) 
Chain 4 Iteration: 7300 / 10000 [ 73%]  (Sampling) 
Chain 4 Iteration: 7400 / 10000 [ 74%]  (Sampling) 
Chain 4 Iteration: 7500 / 10000 [ 75%]  (Sampling) 
Chain 4 Iteration: 7600 / 10000 [ 76%]  (Sampling) 
Chain 4 Iteration: 7700 / 10000 [ 77%]  (Sampling) 
Chain 4 Iteration: 7800 / 10000 [ 78%]  (Sampling) 
Chain 4 Iteration: 7900 / 10000 [ 79%]  (Sampling) 
Chain 4 Iteration: 8000 / 10000 [ 80%]  (Sampling) 
Chain 4 Iteration: 8100 / 10000 [ 81%]  (Sampling) 
Chain 4 Iteration: 8200 / 10000 [ 82%]  (Sampling) 
Chain 4 Iteration: 8300 / 10000 [ 83%]  (Sampling) 
Chain 4 Iteration: 8400 / 10000 [ 84%]  (Sampling) 
Chain 4 Iteration: 8500 / 10000 [ 85%]  (Sampling) 
Chain 4 Iteration: 8600 / 10000 [ 86%]  (Sampling) 
Chain 4 Iteration: 8700 / 10000 [ 87%]  (Sampling) 
Chain 4 Iteration: 8800 / 10000 [ 88%]  (Sampling) 
Chain 4 Iteration: 8900 / 10000 [ 89%]  (Sampling) 
Chain 4 Iteration: 9000 / 10000 [ 90%]  (Sampling) 
Chain 4 Iteration: 9100 / 10000 [ 91%]  (Sampling) 
Chain 4 Iteration: 9200 / 10000 [ 92%]  (Sampling) 
Chain 4 Iteration: 9300 / 10000 [ 93%]  (Sampling) 
Chain 4 Iteration: 9400 / 10000 [ 94%]  (Sampling) 
Chain 4 Iteration: 9500 / 10000 [ 95%]  (Sampling) 
Chain 4 Iteration: 9600 / 10000 [ 96%]  (Sampling) 
Chain 4 Iteration: 9700 / 10000 [ 97%]  (Sampling) 
Chain 4 Iteration: 9800 / 10000 [ 98%]  (Sampling) 
Chain 4 Iteration: 9900 / 10000 [ 99%]  (Sampling) 
Chain 4 Iteration: 10000 / 10000 [100%]  (Sampling) 
Chain 4 finished in 0.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.1 seconds.
Total execution time: 0.6 seconds.

Traceplot

traceplot_ulam(fm)

Rank histogram plot

Shown by Vehtari et al. (2021) to be better than trace plots.

trankplot(fm)

Examining ulam objects

summary(fm)
Inference for Stan model: ulam_cmdstanr_36d52c78ab5e039bffe696e0e0489438-202302200927-1-8b32ac.
4 chains, each with iter=10000; warmup=5000; thin=1; 
post-warmup draws per chain=5000, total post-warmup draws=20000.

       mean se_mean   sd   2.5%   25%   50%   75% 97.5% n_eff Rhat
mu    48.53    0.01 0.71  47.12 48.10 48.53 48.96 49.95  8300    1
sigma  1.91    0.01 0.59   1.12  1.50  1.79  2.18  3.40  7913    1
lp__  -8.22    0.01 1.11 -11.20 -8.66 -7.88 -7.43 -7.13  5780    1

Samples were drawn using NUTS(diag_e) at Mon Feb 20 09:27:21 2023.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

rethinking::precis(): shorter summary

precis(fm)
           mean        sd      5.5%     94.5%    n_eff    Rhat4
mu    48.530352 0.7093896 47.439829 49.634906 8300.361 1.000268
sigma  1.905232 0.5884735  1.206648  2.980274 7913.086 1.000338

n_eff

  • Effect sample size (4 chains \(\times\) 5000 samples = 10,000)
    • We have ~8000
  • Standard deviations are often lower than other parameters
  • Be worried if you have 10,000 samples but n_eff is small
    • Usually accompanied by lots of warnings and high \(\widehat{R}\)

Rhat4

  • Improved version of the Gelman-Rubin \(\widehat{R}\) described by Vehtari et al. (2021)
  • Approaches 1 from above as the chains converge to the same distribution
  • Be concerned if this is more than ~1.01
    • High values (>1.5) are usually accompanied by lots of warnings and low effective sample size.

Sampling gone bad

When sampling fails, it usually fails very badly and very obviously.

  • Investigate your model

Gelman’s Folk Theorem of Statistical Computing:

When you have computational problems, often there’s a problem with your model.

Posteriors

ulam object contains samples:

post <- extract.samples(fm) |> 
  as.data.frame()
head(post)
       mu   sigma
1 47.5197 1.32060
2 50.3552 3.12731
3 47.8664 1.43598
4 48.1501 2.05984
5 48.9498 1.56121
6 47.5633 1.85286
  • Samples are simultaneous at each iteration for each model parameter
    • Important when doing math on samples

Prior vs. Posterior: \(\mu\)

ggplot() +
  geom_line(data = tibble(x = seq(20, 70, length.out = 200),
                          Density = dnorm(x, 45, 10)),
            aes(x, Density),
            color = "#CDDC39",
            linewidth = 2) +
  geom_density(data = post, aes(mu),
               color = "#004C71",
               linewidth = 2) +
  labs(x = expression(mu))

Prior vs. Posterior: \(\mu\)

Prior vs. Posterior: \(\sigma\)

ggplot() +
  geom_line(data = tibble(x = seq(0, 10, length.out = 200),
                          Density = dtruncnorm(x, a = 0, mean = 0, sd = 3)),
            aes(x, Density),
            color = "#CDDC39",
            linewidth = 2) +
  geom_density(data = post, aes(sigma),
               color = "#004C71",
               linewidth = 2) +
  labs(x = expression(sigma))

Prior vs. Posterior: \(\sigma\)

Quantiles

  • Divide the samples into percentiles

89% percentile interval:

quantile(post$mu, probs = c(0.055, 0.945))
    5.5%    94.5% 
47.46445 49.60801

Highest posterior density intervals

  • The highest “mass” region of the distribution

89% HDPI:

HPDI(post$mu)
  |0.89   0.89| 
47.4793 49.6221

Quantiles vs. HPDIs

(Usually nearly) identical for symmetrical distributions:

Quantiles vs. HPDIs

(Can) differ for asymmetrical distributions:

Quantiles vs. HPDIs

  • Use HPDIs if you can.
    • rethinking::HPDI()
    • HDInterval::hdi()
    • coda::HPDinterval()
  • Equivalent for symmetrical (e.g., Gaussian) distributions
  • Better represent the posterior for asymmetrical distributions

References

McElreath, R. 2020. Statistical Rethinking: A Bayesian Course with Examples in R and Stan. 2nd ed. CRC Press.
Vehtari, A., A. Gelman, D. Simpson, B. Carpenter, and P.-C. Bürkner. 2021. Rank-Normalization, Folding, and Localization: An Improved \(\widehat{R}\) for Assessing Convergence of MCMC. Bayesian Analysis 16:667–718.